/*
 * @Descripttion: 全国水土流失代码（带保存），因为归一化太过于费时间，所以暂时注释，但是实际应用必须要归一化，不归一化的话是错的
 * @Author: Haixu He
 * @Date: 2021-12-09 21:17:43
 */
var local = 'Aomen'   //省份
var geometry = ee.FeatureCollection("users/joeyqmf83/Aomen");
Map.addLayer(geometry, {}, "WH_shp")

var dataset = ee.Image("COPERNICUS/Landcover/100m/Proba-V-C3/Global/2019").select('tree-coverfraction').clip(geometry);


//这一部分是处理云盖
function maskS2clouds(image) {
    var qa = image.select('QA60');
    // Bits 10 and 11 are clouds and cirrus, respectively.
    var cloudBitMask = 1 << 10;
    var cirrusBitMask = 1 << 11;

    // Both flags should be set to zero, indicating clear conditions.
    var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
        .and(qa.bitwiseAnd(cirrusBitMask).eq(0));
    return image.updateMask(mask).divide(10000);
}
var dataset_sentinel2 = ee.ImageCollection("COPERNICUS/S2_SR")
    .filterDate('2020-01-01', '2020-12-31')
    .filterBounds(geometry)//
    .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', "less_than", 20)
    .map(maskS2clouds)

var mid_image = dataset_sentinel2.median().clip(geometry);


// // 获取投影信息 然后降采样
// var sentinel2Projection = dataset_sentinel2.first().select('B1').projection();     
// var mid_image = mid_image
//     .reproject({
//       crs: sentinel2Projection
//     })
//     .reduceResolution({
//       reducer: ee.Reducer.mean(),
//       maxPixels: 1024
//     })



var Temp = ee.Image()
//////////////////////////////// Normalization /////////////////////////////////////
function Normalization(image) {
    var minMax = image.reduceRegion({
        reducer: ee.Reducer.minMax(),
        geometry: geometry,
        scale: 30,
        maxPixels: 10e9,
        tileScale: 16
    });
    var temp = ee.ImageCollection.fromImages(
        image.bandNames().map(function (name) {
            name = ee.String(name);
            var band = image.select(name);
            return band.unitScale(ee.Number(minMax.get(name.cat('_min'))), ee.Number(minMax.get(name.cat('_max'))))
        })).toBands().rename(image.bandNames()).clip(geometry);
    return temp
}

///////////////////////////// Yellow ////////////////////////////////////
function Yellow(image) {
    var Y = image.expression(
        "(B3+B4)/2",
        {
            "B4": image.select('B4'),
            "B3": image.select('B3')
        }
    ).rename('Yellow')
    return Y
}
///////////////////////////// NRI ////////////////////////////////////
function NRI(image) {
    var N = image.expression(
        "B8A/B3",
        {
            "B8A": image.select('B8A'),
            "B3": image.select('B3')
        }
    ).rename('NRI')
    return N
}
///////////////////////////// NDSI ////////////////////////////////////
function NDSI(image) {
    var N = image.expression(
        "(B11-B5)/(B11+B5)",
        {
            "B5": image.select('B5'),
            "B11": image.select('B11')
        }
    ).rename('NDSI')

    return N
}

///////////////////////////// NDVI ///////////////////////////////////
var NDVI = mid_image.normalizedDifference(["B8", "B4"]).rename("NDVI")

///////////////////////////// Water Mask ///////////////////////////////////
var mask_water_NDVI = NDVI.updateMask(NDVI.gt(0))//mask water

///////////////////////////// FVC ///////////////////////////////////
function calFVC(BestVI, region, scale) {
    var num = BestVI.reduceRegion({
        reducer: ee.Reducer.percentile([5, 95]),
        geometry: region,
        scale: scale,
        maxPixels: 1e13,
        tileScale: 16
    });
    var min = ee.Number(num.get("NDVI_p5"));
    var max = ee.Number(num.get("NDVI_p95"));
    //quantile and combine
    var greaterPart = BestVI.gt(max);
    var lessPart = BestVI.lt(min);
    var middlePart = ee.Image(1).subtract(greaterPart).subtract(lessPart);
    //calculate FVC
    var tempf1 = BestVI.subtract(min).divide(max.subtract(min));
    var FVC = ee.Image(1).multiply(greaterPart).add(ee.Image(0).multiply(lessPart))
        .add(tempf1.multiply(middlePart))
    return FVC.rename('FVC');
}
var FVC = calFVC(NDVI, geometry, 1000)






///////////////////////////// slope ///////////////////////////////////
var srtm = ee.Image("USGS/SRTMGL1_003");
var dem = ee.Image(srtm).clip(geometry);



// // dem 获取投影信息 然后降采样
// var demProjection = dem.select('elevation').projection();     
// var dem = dem
//     .reproject({
//       crs: demProjection
//     })
//     .reduceResolution({
//       reducer: ee.Reducer.mean(),
//       maxPixels: 900
//     })

var slope = ee.Terrain.slope(dem);


Temp = ee.Image(Yellow(mid_image))
print('Yellow over!', Temp)
Temp = Temp.addBands(NRI(mid_image))
print('NRI over!', Temp)
Temp = Temp.addBands(NDSI(mid_image))
print('NDSI over!', Temp)
Temp = Temp.addBands(NDVI)
print('NDVI over!', Temp)
Temp = Temp.addBands(FVC)
print('FVC over!', Temp)

Temp = Temp.addBands(slope)
print('slope over!', Temp)



// Temp = Normalization(Temp)
// print('Normalization over!',Temp)





////////////////////////////////计算水土流失指数 /////////////////////////////////////
function SEUFM(image) {
    var S = image.expression(
        "(1-FVC)*(1-NRI)*(1-slope)*Yellow*NDSI",
        {
            "FVC": image.select('FVC'),
            "slope": image.select('slope'),
            "NDSI": image.select('NDSI'),
            "Yellow": image.select('Yellow'),
            "NRI": image.select('NRI')
        }
    ).rename('SEUFM')
    return S
}
Temp = Temp.addBands(SEUFM(Temp))
print('SEUFM over!', Temp)
// Temp = Normalization(Temp)
// print('SEUFM Normalization over!',Temp)

// Temp = Temp.where(dataset.lt(30),0)// 筛选
var Temp = Temp.select('SEUFM').updateMask(dataset.gte(20));
print('MASK over!', Temp)


//计算一下各个波段的最大值和最小值
// var minMax = Temp.select('SEUFM').reduceRegion({
//   reducer: ee.Reducer.minMax(),
//   scale: 30,
//   geometry:geometry,
//   maxPixels: 1e13
// });
// print(minMax)



Map.addLayer(Temp.select('SEUFM'), { min: 0, max: 0.6, palette: ["green", "blue", "red"] }, "SEUFM");


//保存tif
Export.image.toDrive({
    image: Temp,
    description: local,
    folder: local,
    scale: 1000,
    region: geometry,
    fileDimensions: 1280,
});

